Florida ITS2 analysis
Nicola G. Kriefall
12/1/21
1 SymPortal pre-submission
Submitting ITS2 samples to Symportal
Prior to submission: retaining only paired reads that match ITS2 primer (we pooled with 16S during sequencing)
# *note: most of this was written by Dr. Carly D. Kenkel
# in Terminal home directory:
# following instructions of installing BBtools from https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/installation-guide/
# 1. download BBMap package, sftp to installation directory
# 2. untar:
tar -xvzf BBMap_(version).tar.gz
# 3. test package:
cd bbmap
~/bin/bbmap/stats.sh in=~/bin/bbmap/resources/phix174_ill.ref.fa.gz
# my adaptors, which I saved as "adaptors.fasta"
# >forward
# AATGATACGGCGACCAC
# >forwardrc
# GTGGTCGCCGTATCATT
# >reverse
# CAAGCAGAAGACGGCATAC
# >reverserc
# GTATGCCGTCTTCTGCTTG
# primers for ITS2:
# >forward
# GTGAATTGCAGAACTCCGTG
# >reverse
# CCTCCGCTTACTTATATGCTT
# making a sample list based on the first phrase before the underscore in the .fastq name
ls *R1_001.fastq | cut -d '_' -f 2 > samples.list
# cuts off the extra words in the .fastq files
for file in $(cat samples.list); do mv ${file}_*R1*.fastq ${file}_R1.fastq; mv ${file}_*R2*.fastq ${file}_R2.fastq; done
# gets rid of reads that still have the adaptor sequence
# [I stopped doing this 'cause it never happened]
#for file in $(cat samples.list); do ~/bin/bbmap/bbduk.sh in1=${file}_R1.fastq in2=${file}_R2.fastq ref=adaptors.fasta out1=${file}_R1_NoIll.fastq out2=${file}_R2_NoIll.fastq; done &>bbduk_NoIll.log
##getting rid of first 4 bases (degenerate primers created them)
for file in $(cat samples.list); do ~/bin/bbmap/bbduk.sh in1=${file}_R1.fastq in2=${file}_R2.fastq ftl=4 out1=${file}_R1_No4N.fastq out2=${file}_R2_No4N.fastq; done &>bbduk_No4N.log
# only keeping reads that start with the primer
for file in $(cat samples.list); do ~/bin/bbmap/bbduk.sh in1=${file}_R1.fastq in2=${file}_R2.fastq k=15 restrictleft=21 literal=GTGAATTGCAGAACTCCGTG,CCTCCGCTTACTTATATGCTT outm1=${file}_R1_ITS.fastq outu1=${file}_R1_check.fastq outm2=${file}_R2_ITS.fastq outu2=${file}_R2_check.fastq; done &>bbduk_ITS.log
# higher k = more reads removed, but can't surpass k=20 or 21
# renamed them to the shorter version again after checking
for file in $(cat samples.list); do mv ${file}_*R1*.fastq ${file}_R1.fastq; mv ${file}_*R2*.fastq ${file}_R2.fastq; done Then transferred all “R1.fastq” & “R2.fastq” files to the folder to be submitted to SymPortal
2 SymPortal ITS2 analysis
Information on Symportal output documents here
Formatting notes:
- Cleaned up file called 188_20211214_03_DBV_20211215T030311.profiles.absolute.abund_and_meta.txt to format for phyloseq
2.1 Setup
## Warning: package 'phyloseq' was built under R version 4.0.3## Loading required package: lattice## Loading required package: plyr## Warning: package 'plyr' was built under R version 4.0.5##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legend## The following object is masked from 'package:plyr':
##
## mutate## Loading required package: permute## Warning: package 'permute' was built under R version 4.0.5## This is vegan 2.6-2## Warning: package 'colorspace' was built under R version 4.0.5##
## Attaching package: 'funfuns'## The following objects are masked from 'package:colorspace':
##
## darken, lighten2.2 Checking for no read samples
counts <- read.csv('flits_symportal_profile_counts.csv',header=TRUE,row.names=1,check.names=FALSE)
#325 samples, 75 type profiles
#+1 sample that's fake 'UK1-O_000' so I can plot UK1-O below
sort(rowSums(counts))## LK1-I_1199_srad_2018 UK2-O_1105_srad_2018 UK1-O_000_srad_2017
## 0 0 1
## UK1-I_1148_ssid_2018 UK1-I_1141_ssid_2018 LK1-I_1192_srad_2018
## 494 550 714
## UK1-O_1173_ssid_2018 LK1-O_1229_srad_2018 UK2-I_1074_srad_2018
## 720 873 1139
## LK1-I_1196_srad_2018 UK2-I_1085_ssid_2018 LK1-I_1193_srad_2018
## 1148 1334 1423
## UK1-I_1144_ssid_2018 UK2-I_1078_srad_2018 UK1-I_1143_ssid_2018
## 1604 1624 1706
## UK2-I_1089_ssid_2018 LK1-I_1197_srad_2018 LK1-I_1191_srad_2018
## 1757 1804 1806
## UK3-I_1037_ssid_2016 UK1-I_1140_ssid_2018 UK2-O_1118_ssid_2018
## 1913 1914 1927
## UK2-I_1084_ssid_2018 UK2-O_1103_srad_2018 LK1-O_1223_srad_2018
## 1958 2005 2039
## LK1-O_1226_srad_2018 LK1-O_1232_ssid_2018 LK1-O_1239_ssid_2018
## 2046 2061 2114
## UK1-O_1176_ssid_2018 UK2-I_1076_srad_2018 UK2-O_1115_ssid_2018
## 2130 2235 2241
## LK1-O_1236_ssid_2018 LK1-I_1194_srad_2018 LK1-O_1221_srad_2018
## 2272 2273 2290
## LK1-O_1233_ssid_2018 UK1-I_1149_ssid_2018 UK3-I_1052_srad_2016
## 2299 2324 2534
## LK1-O_928_srad_2017 UK1-I_1142_ssid_2018 UK3-I_1050_srad_2016
## 2567 2570 2621
## UK2-O_1116_ssid_2018 LK1-I_1190_srad_2018 UK1-O_1174_ssid_2018
## 2685 2755 2764
## UK1-O_1172_ssid_2018 LK1-I_1200_ssid_2018 LK1-O_1235_ssid_2018
## 2822 2830 2893
## LK1-I_1195_srad_2018 LK1-I_1207_ssid_2018 LK1-O_8_ssid_2015
## 2901 2981 3129
## LK1-I_1209_ssid_2018 UK1-I_1147_ssid_2018 UK2-O_97_ssid_2015
## 3134 3223 3306
## UK2-I_1075_srad_2018 LK1-O_1227_srad_2018 UK2-I_1086_ssid_2018
## 3350 3379 3382
## UK2-I_1088_ssid_2018 LK1-I_988_srad_2016 UK2-O_1108_srad_2018
## 3492 3611 3684
## UK1-I_1145_ssid_2018 UK1-O_1171_ssid_2018 UK2-O_1101_srad_2018
## 3815 3876 3907
## UK2-O_1117_ssid_2018 LK1-O_1237_ssid_2018 UK2-O_1119_ssid_2018
## 3911 3918 3995
## UK1-O_1177_ssid_2018 LK1-O_27_srad_2015 LK1-I_1204_ssid_2018
## 4012 4066 4082
## UK2-I_846_srad_2017 UK2-I_1072_srad_2018 LK1-I_990_srad_2016
## 4114 4334 4359
## UK1-O_1170_ssid_2018 UK3-I_1051_srad_2016 LK1-I_1203_ssid_2018
## 4499 4525 4627
## UK2-I_1083_ssid_2018 UK2-O_1107_srad_2018 LK1-O_931_srad_2017
## 4675 4746 4793
## LK1-I_1208_ssid_2018 LK1-O_6_ssid_2015 UK2-I_1071_srad_2018
## 4795 4972 4983
## UK1-I_84_srad_2015 LK1-O_1231_ssid_2018 LK1-I_1202_ssid_2018
## 5039 5058 5148
## LK1-O_1230_ssid_2018 LK1-I_982_ssid_2016 UK2-I_1082_ssid_2018
## 5190 5315 5386
## UK2-I_1081_ssid_2018 LK1-I_1198_srad_2018 UK1-O_783_ssid_2017
## 5409 5417 5432
## UK2-O_1113_ssid_2018 UK2-O_800_ssid_2017 LK1-I_1201_ssid_2018
## 5434 5466 5471
## UK1-I_705_ssid_2017 LK1-I_1206_ssid_2018 UK3-I_1054_srad_2016
## 5487 5731 5836
## UK3-I_1049_srad_2016 LK1-O_1238_ssid_2018 UK2-O_889_ssid_2017
## 5848 5897 5907
## UK1-O_1178_ssid_2018 UK2-O_1109_srad_2018 LK1-O_10_ssid_2015
## 6026 6356 6411
## UK2-O_96_ssid_2015 UK1-O_772_ssid_2017 LK1-O_1222_srad_2018
## 6430 6480 6560
## LK1-O_922_ssid_2017 UK1-O_765_ssid_2017 UK2-O_1112_ssid_2018
## 6589 6621 6844
## LK1-O_26_srad_2015 LK1-O_925_ssid_2017 UK3-O_1026_srad_2016
## 7140 7169 7210
## LK1-O_1228_srad_2018 LK1-O_23_srad_2015 UK1-O_768_ssid_2017
## 7326 7359 7573
## UK2-O_882_ssid_2017 UK2-I_1077_srad_2018 LK1-I_965_ssid_2017
## 7581 7655 7738
## UK2-I_1087_ssid_2018 UK2-I_853_srad_2017 UK2-O_883_ssid_2017
## 7916 7999 8120
## UK2-I_143_srad_2015 LK1-O_24_srad_2015 UK2-I_147_srad_2015
## 8205 8243 8319
## UK2-I_849_srad_2017 UK2-I_852_srad_2017 UK2-O_113_srad_2015
## 8445 8580 8619
## UK2-O_116_srad_2015 LK1-O_30_srad_2015 LK1-O_2_ssid_2015
## 8658 8674 9011
## UK1-O_35_ssid_2015 UK3-I_1039_ssid_2016 LK1-I_981_ssid_2016
## 9045 9081 9090
## UK1-O_1175_ssid_2018 UK1-O_771_ssid_2017 UK2-O_907_srad_2017
## 9100 9127 9153
## LK1-O_4_ssid_2015 UK2-O_1100_srad_2018 UK3-I_1045_ssid_2016
## 9193 9213 9300
## UK2-I_854_srad_2017 LK1-O_3_ssid_2015 UK1-O_763_ssid_2017
## 9571 9707 9798
## UK1-O_769_ssid_2017 UK2-I_850_srad_2017 LK1-I_985_ssid_2016
## 9817 9852 9966
## UK1-I_706_ssid_2017 UK1-I_726_srad_2017 UK2-I_1080_ssid_2018
## 10038 10091 10142
## LK1-I_964_srad_2017 UK1-I_725_ssid_2017 UK2-O_881_ssid_2017
## 10466 10544 10547
## UK3-O_1020_srad_2016 UK2-I_130_ssid_2015 UK2-O_1102_srad_2018
## 10571 10603 10623
## UK2-I_148_srad_2015 UK2-I_825_ssid_2017 UK1-O_762_ssid_2017
## 10632 10812 10875
## UK3-I_1053_srad_2016 UK1-I_707_ssid_2017 UK1-I_731_srad_2017
## 10943 10962 10969
## UK3-I_1048_srad_2016 UK1-O_1179_ssid_2018 LK1-I_978_ssid_2016
## 11057 11267 11647
## UK2-O_92_ssid_2015 UK3-O_1010_ssid_2016 UK1-I_1146_ssid_2018
## 11696 11701 11838
## LK1-I_994_srad_2016 UK2-I_823_ssid_2017 UK1-I_66_ssid_2015
## 11892 12025 12101
## LK1-O_1_ssid_2015 UK1-O_51_srad_2015 LK1-O_929_srad_2017
## 12224 12328 12360
## UK1-I_708_ssid_2017 UK2-O_887_ssid_2017 LK1-I_980_ssid_2016
## 12601 12832 13180
## LK1-O_7_ssid_2015 UK2-I_122_ssid_2015 UK1-I_85_srad_2015
## 13185 13270 13293
## UK3-I_1042_ssid_2016 UK1-I_727_srad_2017 UK2-O_119_srad_2015
## 13311 13317 13392
## UK1-I_701_ssid_2017 UK1-I_733_srad_2017 UK3-I_1040_ssid_2016
## 13398 13429 13466
## LK1-I_986_ssid_2016 UK2-I_851_srad_2017 UK2-O_1106_srad_2018
## 13467 13495 13517
## UK2-I_828_ssid_2017 UK1-O_56_srad_2015 UK2-I_808_ssid_2017
## 13541 13657 13923
## UK1-O_40_ssid_2015 UK2-O_886_ssid_2017 UK2-O_112_srad_2015
## 14077 14106 14293
## LK1-I_949_ssid_2017 UK1-O_33_ssid_2015 UK3-O_1008_ssid_2016
## 14315 14320 14409
## LK1-I_962_srad_2017 LK1-O_943_ssid_2017 UK1-O_38_ssid_2015
## 14467 14502 14652
## UK3-O_1025_srad_2016 UK1-I_732_srad_2017 UK1-I_81_srad_2015
## 14871 15047 15071
## LK1-O_935_srad_2017 UK2-O_885_ssid_2017 UK1-I_704_ssid_2017
## 15076 15391 15520
## UK2-O_909_srad_2017 UK2-I_149_srad_2015 UK3-O_1022_srad_2016
## 15831 15856 15887
## LK1-O_921_ssid_2017 LK1-O_932_srad_2017 LK1-I_977_ssid_2016
## 16000 16277 16476
## LK1-I_993_srad_2016 LK1-I_979_ssid_2016 LK1-O_917_ssid_2017
## 16639 16694 16801
## UK2-O_888_ssid_2017 UK2-O_111_srad_2015 LK1-I_961_srad_2017
## 16802 16832 16915
## UK2-O_120_srad_2015 UK1-O_766_ssid_2017 UK1-I_69_ssid_2015
## 17019 17123 17803
## UK2-O_98_ssid_2015 UK1-O_31_ssid_2015 UK1-O_32_ssid_2015
## 17812 17822 17907
## UK2-I_121_ssid_2015 UK1-I_86_srad_2015 UK1-O_764_ssid_2017
## 17913 18157 18208
## UK1-O_770_ssid_2017 UK1-O_761_ssid_2017 UK2-I_126_ssid_2015
## 18255 18471 18513
## LK1-O_923_ssid_2017 LK1-O_926_srad_2017 UK3-O_1019_srad_2016
## 18540 18626 19010
## LK1-O_22_srad_2015 UK2-O_95_ssid_2015 UK2-I_128_ssid_2015
## 19078 19226 19236
## UK3-O_1018_srad_2016 UK2-I_125_ssid_2015 LK1-I_950_ssid_2017
## 19258 19595 19684
## UK1-O_58_srad_2015 UK2-O_93_ssid_2015 UK1-I_703_ssid_2017
## 19760 20054 20152
## UK2-I_826_ssid_2017 UK2-I_145_srad_2015 UK1-I_709_ssid_2017
## 20261 20332 20403
## LK1-I_948_ssid_2017 UK1-I_65_ssid_2015 LK1-I_984_ssid_2016
## 20408 20800 20846
## UK1-O_36_ssid_2015 LK1-O_919_ssid_2017 UK2-I_824_ssid_2017
## 21174 21265 21277
## UK1-I_68_ssid_2015 LK1-I_951_ssid_2017 UK2-I_123_ssid_2015
## 21408 21520 21552
## UK2-O_915_srad_2017 UK1-I_711_ssid_2017 UK2-O_91_ssid_2015
## 21573 21663 21798
## UK2-O_94_ssid_2015 UK2-I_127_ssid_2015 LK1-O_5_ssid_2015
## 21872 22036 22255
## UK2-I_809_ssid_2017 UK2-O_100_ssid_2015 UK3-I_1041_ssid_2016
## 22499 22568 22603
## LK1-I_963_srad_2017 UK2-I_829_ssid_2017 LK1-O_916_ssid_2017
## 22711 22902 22995
## LK1-O_927_srad_2017 UK3-O_1023_srad_2016 UK1-I_62_ssid_2015
## 23150 23173 23633
## UK1-O_54_srad_2015 UK2-I_129_ssid_2015 LK1-I_976_srad_2017
## 24264 24434 25940
## UK1-I_735_srad_2017 UK2-O_884_ssid_2017 LK1-I_956_srad_2017
## 26183 26322 26615
## UK1-O_37_ssid_2015 UK2-I_124_ssid_2015 UK1-I_82_srad_2015
## 26651 27593 27771
## UK2-O_118_srad_2015 UK1-O_767_ssid_2017 LK1-I_958_srad_2017
## 28102 28250 28449
## UK2-I_855_srad_2017 UK3-I_1046_ssid_2016 UK3-O_1017_srad_2016
## 28604 28623 29057
## UK2-O_913_srad_2017 UK1-O_57_srad_2015 UK3-O_1014_ssid_2016
## 29116 29330 29454
## UK3-O_1011_ssid_2016 LK1-I_947_ssid_2017 UK2-O_890_ssid_2017
## 30004 30367 30503
## UK2-I_807_ssid_2017 UK1-I_64_ssid_2015 UK1-O_39_ssid_2015
## 31352 31518 32161
## UK3-O_1013_ssid_2016 LK1-I_955_ssid_2017 UK2-O_114_srad_2015
## 32219 32224 32661
## UK2-O_798_ssid_2017 UK2-I_827_ssid_2017 UK1-I_63_ssid_2015
## 32936 33128 34025
## UK1-I_702_ssid_2017 UK3-O_1016_ssid_2016 UK2-O_115_srad_2015
## 34417 35721 35763
## LK1-O_21_srad_2015 UK2-O_914_srad_2017 UK1-O_34_ssid_2015
## 36207 37720 38141
## UK2-O_912_srad_2017 UK3-O_1007_ssid_2016 UK2-O_117_srad_2015
## 38907 39109 39575
## UK1-I_70_ssid_2015 UK1-I_61_ssid_2015 UK1-I_67_ssid_2015
## 39625 39692 39935
## UK2-I_821_ssid_2017 UK3-O_1024_srad_2016 UK1-I_710_ssid_2017
## 40254 41155 41230
## UK2-I_822_ssid_2017 LK1-I_983_ssid_2016 UK2-O_908_srad_2017
## 41969 43004 43349
## LK1-I_996_srad_2016 UK3-O_1009_ssid_2016 UK2-I_847_srad_2017
## 44107 44242 44413
## LK1-I_992_srad_2016 UK1-I_89_srad_2015 UK2-O_820_ssid_2017
## 44988 45207 45446
## UK2-I_150_srad_2015 UK3-O_1015_ssid_2016 LK1-I_995_srad_2016
## 45889 46857 47563
## UK2-I_848_srad_2017 UK2-O_99_ssid_2015 UK2-I_830_ssid_2017
## 49830 50651 69634
## UK2-I_144_srad_2015 UK2-I_801_ssid_2017
## 77097 889382.3 Phyloseq object - type profiles
Ran once then saved to re-read in later
# import dataframe holding sample information
samdf <- read.csv("fl_sample_info - ITS2_all.csv")
head(samdf)## sample_num full_site site zone site_zone sample spp year sequencing time
## 1 1 E Sambo LK1 O LK1-O LK1-O_1 ssid 2015 ITS2 Pre
## 2 2 E Sambo LK1 O LK1-O LK1-O_2 ssid 2015 ITS2 Pre
## 3 3 E Sambo LK1 O LK1-O LK1-O_3 ssid 2015 ITS2 Pre
## 4 4 E Sambo LK1 O LK1-O LK1-O_4 ssid 2015 ITS2 Pre
## 5 5 E Sambo LK1 O LK1-O LK1-O_5 ssid 2015 ITS2 Pre
## 6 6 E Sambo LK1 O LK1-O LK1-O_6 ssid 2015 ITS2 Pre
## season
## 1 Spring
## 2 Spring
## 3 Spring
## 4 Spring
## 5 Spring
## 6 Springsamdf$full_sample <- paste0(samdf$sample,"_",samdf$spp,"_",samdf$year)
rownames(samdf) <- samdf$full_sample
# import taxa info
taxa <- read.csv("flits_symportal_taxa.csv",header=TRUE,row.names=1)
#rownames(taxa) <- as.factor(taxa$ITS2_type_profile)
mtaxa <- as.matrix(taxa)
# import counts (absolute abundance from its2 type profiles)
mcounts <- as.matrix(counts.no0)
# Construct phyloseq object
ps <- phyloseq(otu_table(mcounts, taxa_are_rows=FALSE),
sample_data(samdf),
tax_table(mtaxa))
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 75 taxa and 324 samples ]
## sample_data() Sample Data: [ 324 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 75 taxa by 5 taxonomic ranks ]2.5 Bar plots by majority sequence
2.5.1 Assessing type profile distances
Note - used file “20211215T030311_unifrac_profiles_PCoA_coords_B_sqrt.csv”, made the name shorter. I also removed the ‘proportion explained’ line at the end. Also added the first type that shows up as its own column, ‘first_type’
Proportion explained: PC1 0.584828597, PC2 0.146734687 For bray curtis: PC1 0.231, PC2 0.139
library(ggrepel)
setwd("/Volumes/GoogleDrive/My Drive/Flirma/flirma_networks/fl_its2")
uni.cord <- read.csv("unifrac_profiles_PCoA_coords_B_sqrt.csv")
ggplot(uni.cord,aes(x=PC1,y=PC2,color=type_collapsed,label=sample))+
geom_point()+
xlab("PC1 (58.5%)")+
ylab("PC2 (14.7%)")+
#geom_text_repel(max.overlaps=20)+
stat_ellipse()+
theme_bw()## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse## Warning: Removed 3 row(s) containing missing values (geom_path).
#separated B5-1* & B5-2* by the 0 on PC2, except for 1 that was just below 0
#ggsave(file="pca_grouped.pdf",w=12,h=10)
#checking bray curtis
bray.cord <- read.csv("braycurtis_profiles_PCoA_coords_B_sqrt.csv")
ggplot(bray.cord,aes(x=PC1,y=PC2,color=type_bray,label=sample))+
geom_point()+
#geom_text_repel(max.overlaps=20)+
stat_ellipse()## Too few points to calculate an ellipse## Warning: Removed 1 row(s) containing missing values (geom_path).
Important setup
ps.maj <- tax_glom(ps,"Majority_ITS2_collapsed")
#ps.maj <- tax_glom(ps,"Majority_ITS2_sequence_bray")
ntaxa(ps); ntaxa(ps.maj) #18 for unifrac sids, 14 for rads## [1] 75## [1] 14#16 for bray sids
ps.maj <- subset_samples(ps.maj,site!="UK3")
#siderea
ps.maj.sid <- subset_samples(ps.maj,spp=="ssid")
ps.maj.sid.pre <- subset_samples(ps.maj.sid,time=="Pre")
ps.maj.sid.pre.lk1 <- subset_samples(ps.maj.sid.pre,site=="LK1")
ps.maj.sid.pre.uk1 <- subset_samples(ps.maj.sid.pre,site=="UK1")
ps.maj.sid.pre.uk2 <- subset_samples(ps.maj.sid.pre,site=="UK2")
ps.maj.sid.irm <- subset_samples(ps.maj.sid,time=="Irma")
ps.maj.sid.irm.lk1 <- subset_samples(ps.maj.sid.irm,site=="LK1")
ps.maj.sid.irm.uk1 <- subset_samples(ps.maj.sid.irm,site=="UK1")
ps.maj.sid.irm.uk2 <- subset_samples(ps.maj.sid.irm,site=="UK2")
ps.maj.sid.pos <- subset_samples(ps.maj.sid,time=="Post")
ps.maj.sid.pos.lk1 <- subset_samples(ps.maj.sid.pos,site=="LK1")
ps.maj.sid.pos.uk1 <- subset_samples(ps.maj.sid.pos,site=="UK1")
ps.maj.sid.pos.uk2 <- subset_samples(ps.maj.sid.pos,site=="UK2")
#radians
ps.maj.rad <- subset_samples(ps.maj,spp=="srad")
ps.maj.rad.pre <- subset_samples(ps.maj.rad,time=="Pre")
ps.maj.rad.pre.lk1 <- subset_samples(ps.maj.rad.pre,site=="LK1")
ps.maj.rad.pre.uk1 <- subset_samples(ps.maj.rad.pre,site=="UK1")
ps.maj.rad.pre.uk2 <- subset_samples(ps.maj.rad.pre,site=="UK2")
ps.maj.rad.irm <- subset_samples(ps.maj.rad,time=="Irma")
ps.maj.rad.irm.lk1 <- subset_samples(ps.maj.rad.irm,site=="LK1")
ps.maj.rad.irm.uk1 <- subset_samples(ps.maj.rad.irm,site=="UK1")
ps.maj.rad.irm.uk2 <- subset_samples(ps.maj.rad.irm,site=="UK2")
ps.maj.rad.pos <- subset_samples(ps.maj.rad,time=="Post")
ps.maj.rad.pos.lk1 <- subset_samples(ps.maj.rad.pos,site=="LK1")
#ps.maj.rad.pos.uk1 <- subset_samples(ps.maj.rad.pos,site=="UK1")
#not available
ps.maj.rad.pos.uk2 <- subset_samples(ps.maj.rad.pos,site=="UK2")Removing black bars from phyloseq plot
plot_bar2 <- function (physeq, x = "Sample", y = "Abundance", fill = NULL, title = NULL, facet_grid = NULL, border_color = NA)
{
mdf = psmelt(physeq)
p = ggplot(mdf, aes_string(x = x, y = y, fill = fill))
p = p + geom_bar(stat = "identity", position = "stack", color = border_color)
p = p + theme(axis.text.x = element_text(angle = -90, hjust = 0))
if (!is.null(facet_grid)) {
p <- p + facet_grid(facet_grid)
}
if (!is.null(title)) {
p <- p + ggtitle(title)
}
return(p)
}
maj.seq.colors <- c("#474747","#7A7A7A","#A8A8A8","#00294D","#195481","#5882AF","#9AB0CE","#D6DDE9","#D3E0D2","#95B792","#548C4E","#1F5D12","#002F00","#B8773A")2.5.2 The plots
2.5.2.1 S. siderea
LOWER KEYS 1
ps.maj.sid.pre.lk1.rel <- transform_sample_counts(ps.maj.sid.pre.lk1, function(x) x / sum(x))
ps.maj.sid.pre.lk1.rel.zone <- merge_samples(ps.maj.sid.pre.lk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.pre.lk1.zone <- transform_sample_counts(ps.maj.sid.pre.lk1.rel.zone, function(x) x / sum(x))
gg.maj.sid.pre.lk1.zone <- plot_bar2(ps.maj.sid.pre.lk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Pre")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.pre.lk1.zone
ps.maj.sid.irm.lk1.rel <- transform_sample_counts(ps.maj.sid.irm.lk1, function(x) x / sum(x))
ps.maj.sid.irm.lk1.rel.zone <- merge_samples(ps.maj.sid.irm.lk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.irm.lk1.zone <- transform_sample_counts(ps.maj.sid.irm.lk1.rel.zone, function(x) x / sum(x))
gg.maj.sid.irm.lk1.zone <- plot_bar2(ps.maj.sid.irm.lk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Disturb.")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.irm.lk1.zone
ps.maj.sid.pos.lk1.rel <- transform_sample_counts(ps.maj.sid.pos.lk1, function(x) x / sum(x))
ps.maj.sid.pos.lk1.rel.zone <- merge_samples(ps.maj.sid.pos.lk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.pos.lk1.zone <- transform_sample_counts(ps.maj.sid.pos.lk1.rel.zone, function(x) x / sum(x))
gg.maj.sid.pos.lk1.zone <- plot_bar2(ps.maj.sid.pos.lk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Post")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.pos.lk1.zone
UPPER KEYS 1
ps.maj.sid.pre.uk1.rel <- transform_sample_counts(ps.maj.sid.pre.uk1, function(x) x / sum(x))
ps.maj.sid.pre.uk1.rel.zone <- merge_samples(ps.maj.sid.pre.uk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.pre.uk1.zone <- transform_sample_counts(ps.maj.sid.pre.uk1.rel.zone, function(x) x / sum(x))
gg.maj.sid.pre.uk1.zone <- plot_bar2(ps.maj.sid.pre.uk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Pre")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.pre.uk1.zone
ps.maj.sid.irm.uk1.rel <- transform_sample_counts(ps.maj.sid.irm.uk1, function(x) x / sum(x))
ps.maj.sid.irm.uk1.rel.zone <- merge_samples(ps.maj.sid.irm.uk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.irm.uk1.zone <- transform_sample_counts(ps.maj.sid.irm.uk1.rel.zone, function(x) x / sum(x))
gg.maj.sid.irm.uk1.zone <- plot_bar2(ps.maj.sid.irm.uk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Disturb.")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.irm.uk1.zone
ps.maj.sid.pos.uk1.rel <- transform_sample_counts(ps.maj.sid.pos.uk1, function(x) x / sum(x))
ps.maj.sid.pos.uk1.rel.zone <- merge_samples(ps.maj.sid.pos.uk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.pos.uk1.zone <- transform_sample_counts(ps.maj.sid.pos.uk1.rel.zone, function(x) x / sum(x))
gg.maj.sid.pos.uk1.zone <- plot_bar2(ps.maj.sid.pos.uk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Post")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.pos.uk1.zone
UPPER KEYS 2
ps.maj.sid.pre.uk2.rel <- transform_sample_counts(ps.maj.sid.pre.uk2, function(x) x / sum(x))
ps.maj.sid.pre.uk2.rel.zone <- merge_samples(ps.maj.sid.pre.uk2.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.pre.uk2.zone <- transform_sample_counts(ps.maj.sid.pre.uk2.rel.zone, function(x) x / sum(x))
gg.maj.sid.pre.uk2.zone <- plot_bar2(ps.maj.sid.pre.uk2.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Pre")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.pre.uk2.zone
ps.maj.sid.irm.uk2.rel <- transform_sample_counts(ps.maj.sid.irm.uk2, function(x) x / sum(x))
ps.maj.sid.irm.uk2.rel.zone <- merge_samples(ps.maj.sid.irm.uk2.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.irm.uk2.zone <- transform_sample_counts(ps.maj.sid.irm.uk2.rel.zone, function(x) x / sum(x))
gg.maj.sid.irm.uk2.zone <- plot_bar2(ps.maj.sid.irm.uk2.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Disturb.")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.irm.uk2.zone
ps.maj.sid.pos.uk2.rel <- transform_sample_counts(ps.maj.sid.pos.uk2, function(x) x / sum(x))
ps.maj.sid.pos.uk2.rel.zone <- merge_samples(ps.maj.sid.pos.uk2.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.pos.uk2.zone <- transform_sample_counts(ps.maj.sid.pos.uk2.rel.zone, function(x) x / sum(x))
gg.maj.sid.pos.uk2.zone <- plot_bar2(ps.maj.sid.pos.uk2.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Post")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.sid.pos.uk2.zone
All the sites
ggarrange(gg.maj.sid.pre.lk1.zone,gg.maj.sid.irm.lk1.zone,gg.maj.sid.pos.lk1.zone,gg.maj.sid.pre.uk1.zone,gg.maj.sid.irm.uk1.zone,gg.maj.sid.pos.uk1.zone,gg.maj.sid.pre.uk2.zone,gg.maj.sid.irm.uk2.zone,gg.maj.sid.pos.uk2.zone,nrow=3,ncol=3,common.legend=T,legend="right")
2.5.2.2 S. radians
LOWER KEYS 1
ps.maj.rad.pre.lk1.rel <- transform_sample_counts(ps.maj.rad.pre.lk1, function(x) x / sum(x))
ps.maj.rad.pre.lk1.rel.zone <- merge_samples(ps.maj.rad.pre.lk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.rad.pre.lk1.zone <- transform_sample_counts(ps.maj.rad.pre.lk1.rel.zone, function(x) x / sum(x))
gg.maj.rad.pre.lk1.zone <- plot_bar2(ps.maj.rad.pre.lk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Pre")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.rad.pre.lk1.zone
ps.maj.rad.irm.lk1.rel <- transform_sample_counts(ps.maj.rad.irm.lk1, function(x) x / sum(x))
ps.maj.rad.irm.lk1.rel.zone <- merge_samples(ps.maj.rad.irm.lk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.rad.irm.lk1.zone <- transform_sample_counts(ps.maj.rad.irm.lk1.rel.zone, function(x) x / sum(x))
gg.maj.rad.irm.lk1.zone <- plot_bar2(ps.maj.rad.irm.lk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Disturb.")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.rad.irm.lk1.zone
ps.maj.rad.pos.lk1.rel <- transform_sample_counts(ps.maj.rad.pos.lk1, function(x) x / sum(x))
ps.maj.rad.pos.lk1.rel.zone <- merge_samples(ps.maj.rad.pos.lk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.rad.pos.lk1.zone <- transform_sample_counts(ps.maj.rad.pos.lk1.rel.zone, function(x) x / sum(x))
gg.maj.rad.pos.lk1.zone <- plot_bar2(ps.maj.rad.pos.lk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Post")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.rad.pos.lk1.zone
UPPER KEYS 1
ps.maj.rad.pre.uk1.rel <- transform_sample_counts(ps.maj.rad.pre.uk1, function(x) x / sum(x))
ps.maj.rad.pre.uk1.rel.zone <- merge_samples(ps.maj.rad.pre.uk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.rad.pre.uk1.zone <- transform_sample_counts(ps.maj.rad.pre.uk1.rel.zone, function(x) x / sum(x))
gg.maj.rad.pre.uk1.zone <- plot_bar2(ps.maj.rad.pre.uk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Pre")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.rad.pre.uk1.zone
ps.maj.rad.irm.uk1.rel <- transform_sample_counts(ps.maj.rad.irm.uk1, function(x) x / sum(x))
ps.maj.rad.irm.uk1.rel.zone <- merge_samples(ps.maj.rad.irm.uk1.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.rad.irm.uk1.zone <- transform_sample_counts(ps.maj.rad.irm.uk1.rel.zone, function(x) x / sum(x))
gg.maj.rad.irm.uk1.zone <- plot_bar2(ps.maj.rad.irm.uk1.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Disturb.")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.rad.irm.uk1.zone
#uk1 not available for Post
#ggarrange(gg.maj.rad.pre.uk1.zone,gg.maj.rad.irm.uk1.zone,gg.maj.rad.pos.uk1.zone,nrow=1)UPPER KEYS 2
ps.maj.rad.pre.uk2.rel <- transform_sample_counts(ps.maj.rad.pre.uk2, function(x) x / sum(x))
ps.maj.rad.pre.uk2.rel.zone <- merge_samples(ps.maj.rad.pre.uk2.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.rad.pre.uk2.zone <- transform_sample_counts(ps.maj.rad.pre.uk2.rel.zone, function(x) x / sum(x))
gg.maj.rad.pre.uk2.zone <- plot_bar2(ps.maj.rad.pre.uk2.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Pre")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.rad.pre.uk2.zone
ps.maj.rad.irm.uk2.rel <- transform_sample_counts(ps.maj.rad.irm.uk2, function(x) x / sum(x))
ps.maj.rad.irm.uk2.rel.zone <- merge_samples(ps.maj.rad.irm.uk2.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.rad.irm.uk2.zone <- transform_sample_counts(ps.maj.rad.irm.uk2.rel.zone, function(x) x / sum(x))
gg.maj.rad.irm.uk2.zone <- plot_bar2(ps.maj.rad.irm.uk2.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Disturb.")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.rad.irm.uk2.zone
ps.maj.rad.pos.uk2.rel <- transform_sample_counts(ps.maj.rad.pos.uk2, function(x) x / sum(x))
ps.maj.rad.pos.uk2.rel.zone <- merge_samples(ps.maj.rad.pos.uk2.rel, "zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.rad.pos.uk2.zone <- transform_sample_counts(ps.maj.rad.pos.uk2.rel.zone, function(x) x / sum(x))
gg.maj.rad.pos.uk2.zone <- plot_bar2(ps.maj.rad.pos.uk2.zone, fill="Majority_ITS2_collapsed")+
ggtitle("Post")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_x_discrete(labels=c("Inshore","Offshore"))+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
theme(axis.text.x=element_text(angle=45))
gg.maj.rad.pos.uk2.zone
All the sites
ggarrange(gg.maj.rad.pre.lk1.zone,gg.maj.rad.irm.lk1.zone,gg.maj.rad.pos.lk1.zone,gg.maj.rad.pre.uk1.zone,gg.maj.rad.irm.uk1.zone,gg.maj.rad.pre.uk2.zone,gg.maj.rad.irm.uk2.zone,gg.maj.rad.pos.uk2.zone,nrow=3,ncol=3,common.legend=T,legend="right")
2.5.3 Different version of the same type profile plot
ps.maj.sid.pre.rel <- transform_sample_counts(ps.maj.sid.pre, function(x) x / sum(x))
ps.maj.sid.pre.rel.sz <- merge_samples(ps.maj.sid.pre.rel, "site_zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.pre.rel.sz <- transform_sample_counts(ps.maj.sid.pre.rel.sz, function(x) x / sum(x))
gg.bar.pre <- plot_bar2(ps.maj.sid.pre.rel.sz, fill="Majority_ITS2_collapsed")+
ggtitle("Pre")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_fill_manual(values=maj.seq.colors)+
geom_bar(stat = "identity")
gg.bar.pre
ps.maj.sid.irm.rel <- transform_sample_counts(ps.maj.sid.irm, function(x) x / sum(x))
ps.maj.sid.irm.rel.sz <- merge_samples(ps.maj.sid.irm.rel, "site_zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.irm.rel.sz <- transform_sample_counts(ps.maj.sid.irm.rel.sz, function(x) x / sum(x))
gg.bar.irm <- plot_bar(ps.maj.sid.irm.rel.sz, fill="Majority_ITS2_collapsed")+
ggtitle("Disturbance")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_fill_manual(values=maj.seq.colors)
gg.bar.irm
ps.maj.sid.pos.rel <- transform_sample_counts(ps.maj.sid.pos, function(x) x / sum(x))
ps.maj.sid.pos.rel.sz <- merge_samples(ps.maj.sid.pos.rel, "site_zone")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionps.maj.sid.pos.rel.sz <- transform_sample_counts(ps.maj.sid.pos.rel.sz, function(x) x / sum(x))
gg.bar.pos <- plot_bar(ps.maj.sid.pos.rel.sz, fill="Majority_ITS2_collapsed")+
ggtitle("Post")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_fill_manual(values=maj.seq.colors)
gg.bar.pos
gg.bar.its2 <- ggarrange(gg.bar.pre,gg.bar.irm,gg.bar.pos,ncol=3,common.legend=T,legend="right")
gg.bar.its2
Type profile - Bray-Curtis distance
maj.seq.colors <- c("#1B1B1B","#969696","#F9F9F9","#00294D","#195481","#5882AF","#9AB0CE","#D6DDE9","#D3E0D2","#95B792","#548C4E","#1F5D12","#002F00","#B8773A")
ps.maj.sid.pre.rel <- transform_sample_counts(ps.maj.sid.pre, function(x) x / sum(x))
ps.maj.sid.pre.rel.sz <- merge_samples(ps.maj.sid.pre.rel, "site_zone")
ps.maj.sid.pre.rel.sz <- transform_sample_counts(ps.maj.sid.pre.rel.sz, function(x) x / sum(x))
gg.bar.pre <- plot_bar(ps.maj.sid.pre.rel.sz, fill="Majority_ITS2_sequence_bray")+
ggtitle("Pre")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_fill_manual(values=maj.seq.colors)
gg.bar.pre
ps.maj.sid.irm.rel <- transform_sample_counts(ps.maj.sid.irm, function(x) x / sum(x))
ps.maj.sid.irm.rel.sz <- merge_samples(ps.maj.sid.irm.rel, "site_zone")
ps.maj.sid.irm.rel.sz <- transform_sample_counts(ps.maj.sid.irm.rel.sz, function(x) x / sum(x))
gg.bar.irm <- plot_bar(ps.maj.sid.irm.rel.sz, fill="Majority_ITS2_sequence_bray")+
ggtitle("Disturbance")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_fill_manual(values=maj.seq.colors)
gg.bar.irm
ps.maj.sid.pos.rel <- transform_sample_counts(ps.maj.sid.pos, function(x) x / sum(x))
ps.maj.sid.pos.rel.sz <- merge_samples(ps.maj.sid.pos.rel, "site_zone")
ps.maj.sid.pos.rel.sz <- transform_sample_counts(ps.maj.sid.pos.rel.sz, function(x) x / sum(x))
gg.bar.pos <- plot_bar(ps.maj.sid.pos.rel.sz, fill="Majority_ITS2_sequence_bray")+
ggtitle("Post")+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
xlab("")+
scale_fill_manual(values=maj.seq.colors)
gg.bar.pos
gg.bar.its2 <- ggarrange(gg.bar.pre,gg.bar.irm,gg.bar.pos,ncol=3,common.legend=T,legend="right")
gg.bar.its2
#ggsave("bar.its2.sids.pdf",width=9,height=4)2.6 Bar plot by type profile
2.6.1 Setup - packages
#remotes::install_github("KarstensLab/microshades")
library("microshades")
#remotes::install_github("mikemc/speedyseq")
library("speedyseq")##
## Attaching package: 'speedyseq'## The following objects are masked from 'package:phyloseq':
##
## filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom,
## tip_glom, transform_sample_counts2.6.2 Setup - data
2.6.3 Microshades!
2.6.3.1 Microshades plot - sids
# Use microshades function prep_mdf to agglomerate, normalize, and melt the phyloseq object
mdf.ps.sid <- prep_mdf(ps.sid.no0,subgroup_level="Majority_ITS2_collapsed")
type_groups <- c(unique(mdf.ps.sid$Clade))
# Create a color object for the specified data
col.mdf.ps.sid <- create_color_dfs(mdf.ps.sid, top_orientation = FALSE,group_level="Clade",subgroup_level="Majority_ITS2_collapsed",selected_groups=c("A","B","C","D"),cvd=TRUE)
# Extract
col.mdf.ps.sid.m <- col.mdf.ps.sid$mdf
col.mdf.ps.sid.c <- col.mdf.ps.sid$cdf
col.mdf.ps.sid.m$zone <- sub("I","Inshore",col.mdf.ps.sid.m$zone)
col.mdf.ps.sid.m$zone <- sub("O","Offshore",col.mdf.ps.sid.m$zone)
col.mdf.ps.sid.m$time <- factor(col.mdf.ps.sid.m$time,levels=c("Pre","Irma","Post"))
col.mdf.ps.sid.c.new <- color_reassign(col.mdf.ps.sid.c,
group_assignment = c("A","B","C","D"),
color_assignment = c("micro_gray","micro_cvd_blue","micro_cvd_green","micro_brown"),group_level="Clade")
col.mdf.ps.sid.m.lk1 <- subset(col.mdf.ps.sid.m,site=="LK1")
ms.plot.lk1 <- plot_microshades(col.mdf.ps.sid.m.lk1, col.mdf.ps.sid.c.new)+
facet_wrap(zone~time,scales="free_x")+
theme_cowplot()+
theme(axis.text.x=element_blank())+
ggtitle("Lower Keys 1")+
guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence"))
#theme(axis.text.x=element_text(angle=45,hjust=1))
ms.plot.lk1
col.mdf.ps.sid.m.uk1 <- subset(col.mdf.ps.sid.m,site=="UK1")
ms.plot.uk1 <- plot_microshades(col.mdf.ps.sid.m.uk1, col.mdf.ps.sid.c.new)+
facet_wrap(zone~time,scales="free_x")+
theme_cowplot()+
theme(axis.text.x=element_blank())+
ggtitle("Upper Keys 1")+
guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence"))
ms.plot.uk1
col.mdf.ps.sid.m.uk2 <- subset(col.mdf.ps.sid.m,site=="UK2")
ms.plot.uk2 <- plot_microshades(col.mdf.ps.sid.m.uk2, col.mdf.ps.sid.c.new)+
facet_wrap(zone~time,scales="free_x")+
theme_cowplot()+
theme(axis.text.x=element_blank())+
ggtitle("Upper Keys 2")+
guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence"))
ms.plot.uk2
ms.plot.sid <- ggarrange(ms.plot.lk1,ms.plot.uk1,ms.plot.uk2,nrow=3,ncol=1,labels="AUTO")
annotate_figure(ms.plot.sid,top=text_grob("S. siderea",face="italic",size=16))
2.6.3.2 Microshades plot - rads
# Use microshades function prep_mdf to agglomerate, normalize, and melt the phyloseq object
mdf.ps.rad <- prep_mdf(ps.rad.no0,subgroup_level="Majority_ITS2_collapsed")
type_groups <- c(unique(mdf.ps.rad$Clade))
# Create a color object for the specified data
col.mdf.ps.rad <- create_color_dfs(mdf.ps.rad, top_orientation = FALSE,group_level="Clade",subgroup_level="Majority_ITS2_collapsed",selected_groups=c("A","B","C","D"),cvd=TRUE)
# Extract
col.mdf.ps.rad.m <- col.mdf.ps.rad$mdf
col.mdf.ps.rad.c <- col.mdf.ps.rad$cdf
col.mdf.ps.rad.m$zone <- sub("I","Inshore",col.mdf.ps.rad.m$zone)
col.mdf.ps.rad.m$zone <- sub("O","Offshore",col.mdf.ps.rad.m$zone)
col.mdf.ps.rad.m$time <- factor(col.mdf.ps.rad.m$time,levels=c("Pre","Irma","Post"))
col.mdf.ps.rad.c.new <- color_reassign(col.mdf.ps.rad.c,
group_assignment = c("A","B","C","D"),
color_assignment = c("micro_gray","micro_cvd_blue","micro_cvd_green","micro_brown"),group_level="Clade")
col.mdf.ps.rad.m.lk1 <- subset(col.mdf.ps.rad.m,site=="LK1")
ms.plot.lk1 <- plot_microshades(col.mdf.ps.rad.m.lk1, col.mdf.ps.rad.c.new)+
facet_wrap(zone~time,scales="free_x")+
theme_cowplot()+
theme(axis.text.x=element_blank())+
ggtitle("Lower Keys 1")+
guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence"))
ms.plot.lk1
col.mdf.ps.rad.m.uk1 <- subset(col.mdf.ps.rad.m,site=="UK1")
ms.plot.uk1 <- plot_microshades(col.mdf.ps.rad.m.uk1, col.mdf.ps.rad.c.new)+
facet_wrap(zone~time,scales="free_x")+
theme_cowplot()+
theme(axis.text.x=element_blank())+
ggtitle("Upper Keys 1")+
guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence"))
ms.plot.uk1
col.mdf.ps.rad.m.uk2 <- subset(col.mdf.ps.rad.m,site=="UK2")
ms.plot.uk2 <- plot_microshades(col.mdf.ps.rad.m.uk2, col.mdf.ps.rad.c.new)+
facet_wrap(zone~time,scales="free_x")+
theme_cowplot()+
theme(axis.text.x=element_blank())+
ggtitle("Upper Keys 1")+
guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence"))
ms.plot.uk2
ms.plot.rad <- ggarrange(ms.plot.lk1,ms.plot.uk1,ms.plot.uk2,nrow=3,ncol=1)
annotate_figure(ms.plot.rad,top=text_grob("S. radians",face="italic",size=16))
3 PCoAs/stats
3.1 Setup
library(cowplot)
library(phyloseq)
library(vegan)
library(ggplot2)
setwd("/Volumes/GoogleDrive/My Drive/Flirma/flirma_networks/fl_its2")
ps.pm <- readRDS("ps.its2.postmed.rds")Run once, then skip
counts <- read.csv('flits_postmed.csv',header=TRUE,row.names=1,check.names=FALSE)
#325 samples, 1268 post med seqs
#in relative abundance format
# import dataframe holding sample information
samdf <- read.csv("fl_sample_info - ITS2_all.csv")
head(samdf)
samdf$full_sample <- paste0(samdf$sample,"_",samdf$spp,"_",samdf$year)
rownames(samdf) <- samdf$full_sample
# import counts (relative abundance from its2 type profiles)
mcounts <- as.matrix(counts)
# Construct phyloseq object
ps.pm1 <- phyloseq(otu_table(mcounts, taxa_are_rows=FALSE),
sample_data(samdf))
ps.pm1
ps.pm <- subset_samples(ps.pm1,site!="UK3")
ps.pm
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 1268 taxa and 293 samples ]
# sample_data() Sample Data: [ 293 samples by 12 sample variables ]
#saveRDS(ps.pm,file="ps.its2.postmed.rds")3.2 Stats by time - sids
3.2.1 Setup
ps.pm.sid <- subset_samples(ps.pm,spp=="ssid")
##for time stats
ps.pm.sid.lk1 <- subset_samples(ps.pm.sid,site=="LK1")
ps.pm.sid.uk1 <- subset_samples(ps.pm.sid,site=="UK1")
ps.pm.sid.uk2 <- subset_samples(ps.pm.sid,site=="UK2")
ps.pm.sid.lk1i <- subset_samples(ps.pm.sid.lk1,zone=="I")
ps.pm.sid.lk1o <- subset_samples(ps.pm.sid.lk1,zone=="O")
ps.pm.sid.uk1i <- subset_samples(ps.pm.sid.uk1,zone=="I")
ps.pm.sid.uk1o <- subset_samples(ps.pm.sid.uk1,zone=="O")
ps.pm.sid.uk2i <- subset_samples(ps.pm.sid.uk2,zone=="I")
ps.pm.sid.uk2o <- subset_samples(ps.pm.sid.uk2,zone=="O")3.2.2 Time - LK1 sids
##time - inshore
seq.sid.lk1i <- data.frame(ps.pm.sid.lk1i@otu_table)
samdf.sid.lk1i <- data.frame(ps.pm.sid.lk1i@sam_data)
dist.sid.lk1i <- vegdist(seq.sid.lk1i)
bet.sid.lk1i <- betadisper(dist.sid.lk1i,samdf.sid.lk1i$time)
#anova(bet.ps.pm)
permutest(bet.sid.lk1i,pairwise=TRUE,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.08061 0.040304 0.9945 999 0.375
## Residuals 23 0.93215 0.040528
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.029000 0.590
## Post 0.034299 0.446
## Pre 0.615872 0.442412#post & Irma differ
# Irma Post Pre
# Irma 0.032000 0.603
# Post 0.034299 0.427
# Pre 0.615872 0.442412
plot(bet.sid.lk1i) 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.sid.lk1i ~ time, data = samdf.sid.lk1i, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 1.8496 0.25275 3.8898 0.004 **
## Residual 23 5.4683 0.74725
## Total 25 7.3179 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Post vs Irma 3.6451374 0.20658028 0.002 0.002
## 2 Post vs Pre 7.7147744 0.31215233 0.001 0.001
## 3 Irma vs Pre 0.7382707 0.04690926 0.441 0.441##post different than other two
# pairs F.Model R2 p.value p.adjusted
# 1 Post vs Irma 3.6451374 0.20658028 0.004 0.004
# 2 Post vs Pre 7.7147744 0.31215233 0.002 0.002
# 3 Irma vs Pre 0.7382707 0.04690926 0.430 0.430
##time - offshore
seq.sid.lk1o <- data.frame(ps.pm.sid.lk1o@otu_table)
samdf.sid.lk1o <- data.frame(ps.pm.sid.lk1o@sam_data)
dist.sid.lk1o <- vegdist(seq.sid.lk1o)
bet.sid.lk1o <- betadisper(dist.sid.lk1o,samdf.sid.lk1o$time)
#anova(bet.ps.pm)
permutest(bet.sid.lk1o,pairwise=TRUE,permutations=999)#ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.12883 0.064413 0.8942 999 0.433
## Residuals 23 1.65672 0.072031
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.86900 0.310
## Post 0.88312 0.341
## Pre 0.31083 0.33073
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.sid.lk1o ~ time, data = samdf.sid.lk1o, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 2.0304 0.27913 4.4531 0.006 **
## Residual 23 5.2434 0.72087
## Total 25 7.2738 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Post vs Pre 7.621797 0.3226595 0.001 0.001
## 2 Post vs Irma 3.638140 0.1951987 0.021 0.021
## 3 Pre vs Irma 2.032004 0.1193050 0.139 0.1393.2.3 Time - UK1 sids
##time
seq.sid.uk1i <- data.frame(ps.pm.sid.uk1i@otu_table)
samdf.sid.uk1i <- data.frame(ps.pm.sid.uk1i@sam_data)
dist.sid.uk1i <- vegdist(seq.sid.uk1i)
bet.sid.uk1i <- betadisper(dist.sid.uk1i,samdf.sid.uk1i$time)
#anova(bet.ps.pm)
permutest(bet.sid.uk1i,pairwise=TRUE,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.0743 0.037151 0.9067 999 0.434
## Residuals 29 1.1883 0.040975
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.33200 0.684
## Post 0.30524 0.232
## Pre 0.69047 0.20993
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.sid.uk1i ~ time, data = samdf.sid.uk1i, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 0.6774 0.12238 2.022 0.041 *
## Residual 29 4.8576 0.87762
## Total 31 5.5349 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Post vs Pre 2.804997 0.13482324 0.007 0.007
## 2 Post vs Irma 2.497810 0.11102459 0.037 0.037
## 3 Pre vs Irma 0.604099 0.02931936 0.568 0.568# pairs F.Model R2 p.value p.adjusted
# 1 Post vs Pre 2.804997 0.13482324 0.011 0.011
# 2 Post vs Irma 2.497810 0.11102459 0.030 0.030
# 3 Pre vs Irma 0.604099 0.02931936 0.537 0.537
#post different from other two
##time - offshore
seq.sid.uk1o <- data.frame(ps.pm.sid.uk1o@otu_table)
samdf.sid.uk1o <- data.frame(ps.pm.sid.uk1o@sam_data)
dist.sid.uk1o <- vegdist(seq.sid.uk1o)
bet.sid.uk1o <- betadisper(dist.sid.uk1o,samdf.sid.uk1o$time)
#anova(bet.ps.pm)
permutest(bet.sid.uk1o,pairwise=TRUE,permutations=999)#ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.05342 0.026708 0.3119 999 0.742
## Residuals 30 2.56879 0.085626
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.96000 0.525
## Post 0.96642 0.501
## Pre 0.51116 0.50881
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.sid.uk1o ~ time, data = samdf.sid.uk1o, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 0.7656 0.10339 1.7297 0.125
## Residual 30 6.6393 0.89661
## Total 32 7.4049 1.00000## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Post vs Pre 2.3684810 0.11628167 0.076 0.076
## 2 Post vs Irma 2.3930846 0.10229881 0.081 0.081
## 3 Pre vs Irma 0.5068001 0.02356465 0.465 0.4653.2.4 Time - UK2 sids
##time - inshore
seq.sid.uk2i <- data.frame(ps.pm.sid.uk2i@otu_table)
samdf.sid.uk2i <- data.frame(ps.pm.sid.uk2i@sam_data)
dist.sid.uk2i <- vegdist(seq.sid.uk2i)
bet.sid.uk2i <- betadisper(dist.sid.uk2i,samdf.sid.uk2i$time)
#anova(bet.ps.pm)
permutest(bet.sid.uk2i,pairwise=TRUE,permutations=999)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.21882 0.109412 5.6276 999 0.009 **
## Residuals 31 0.60270 0.019442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.055000 0.049
## Post 0.057523 0.015
## Pre 0.050815 0.013228##Pre & irma are different
# Irma Post Pre
# Irma 0.439000 0.016
# Post 0.431974 0.187
# Pre 0.017649 0.188058
plot(bet.sid.uk2i) #pre less constrained
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.sid.uk2i ~ time, data = samdf.sid.uk2i, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 1.1924 0.17282 3.2385 0.001 ***
## Residual 31 5.7073 0.82718
## Total 33 6.8997 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Post vs Irma 3.494325 0.13706285 0.001 0.001
## 2 Post vs Pre 3.981944 0.18114612 0.002 0.002
## 3 Irma vs Pre 1.932622 0.08075262 0.120 0.120##post different than other two
# pairs F.Model R2 p.value p.adjusted
# 1 Pre vs Post 1.958299 0.05766776 0.028 0.028
# 2 Pre vs Irma 1.398358 0.03461423 0.164 0.164
# 3 Post vs Irma 1.613048 0.03971747 0.050 0.050
##time - offshore
seq.sid.uk2o <- data.frame(ps.pm.sid.uk2o@otu_table)
samdf.sid.uk2o <- data.frame(ps.pm.sid.uk2o@sam_data)
dist.sid.uk2o <- vegdist(seq.sid.uk2o)
bet.sid.uk2o <- betadisper(dist.sid.uk2o,samdf.sid.uk2o$time)
#anova(bet.ps.pm)
permutest(bet.sid.uk2o,pairwise=TRUE,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.08429 0.042143 0.8818 999 0.428
## Residuals 27 1.29039 0.047792
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.31900 0.804
## Post 0.33697 0.110
## Pre 0.77312 0.11717
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.sid.uk2o ~ time, data = samdf.sid.uk2o, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 0.6775 0.09085 1.349 0.226
## Residual 27 6.7806 0.90915
## Total 29 7.4581 1.00000## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Post vs Irma 1.6446957 0.08372212 0.116 0.116
## 2 Post vs Pre 2.0931451 0.12245524 0.095 0.095
## 3 Irma vs Pre 0.6640866 0.03065380 0.471 0.4713.3 Plot by time - sids
3.3.1 PCOA LK1 sids
ps.pm.sid.lk1i.no0 <- prune_taxa(taxa_sums(ps.pm.sid.lk1i)!=0,ps.pm.sid.lk1i)
ord.pm.sid.lk1i.no0 <- ordinate(ps.pm.sid.lk1i.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.sid.lk1i.no0, ord.pm.sid.lk1i.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
ps.pm.sid.lk1o.no0 <- prune_taxa(taxa_sums(ps.pm.sid.lk1o)!=0,ps.pm.sid.lk1o)
ord.pm.sid.lk1o.no0 <- ordinate(ps.pm.sid.lk1o.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.sid.lk1o.no0, ord.pm.sid.lk1o.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
3.3.2 PCOA UK1 sids
ps.pm.sid.uk1i.no0 <- prune_taxa(taxa_sums(ps.pm.sid.uk1i)!=0,ps.pm.sid.uk1i)
ord.pm.sid.uk1i.no0 <- ordinate(ps.pm.sid.uk1i.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.sid.uk1i.no0, ord.pm.sid.uk1i.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()
ps.pm.sid.uk1o.no0 <- prune_taxa(taxa_sums(ps.pm.sid.uk1o)!=0,ps.pm.sid.uk1o)
ord.pm.sid.uk1o.no0 <- ordinate(ps.pm.sid.uk1o.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.sid.uk1o.no0, ord.pm.sid.uk1o.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
3.3.3 PCOA UK2 sids
ps.pm.sid.uk2i.no0 <- prune_taxa(taxa_sums(ps.pm.sid.uk2i)!=0,ps.pm.sid.uk2i)
ord.pm.sid.uk2i.no0 <- ordinate(ps.pm.sid.uk2i.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.sid.uk2i.no0, ord.pm.sid.uk2i.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()
ps.pm.sid.uk2o.no0 <- prune_taxa(taxa_sums(ps.pm.sid.uk2o)!=0,ps.pm.sid.uk2o)
ord.pm.sid.uk2o.no0 <- ordinate(ps.pm.sid.uk2o.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.sid.uk2o.no0, ord.pm.sid.uk2o.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()
3.4 Stats by time - rads
3.4.1 Setup
ps.pm.rad <- subset_samples(ps.pm,spp=="srad")
##for time stats
ps.pm.rad.lk1 <- subset_samples(ps.pm.rad,site=="LK1")
ps.pm.rad.uk1 <- subset_samples(ps.pm.rad,site=="UK1")
ps.pm.rad.uk2 <- subset_samples(ps.pm.rad,site=="UK2")
ps.pm.rad.lk1i <- subset_samples(ps.pm.rad.lk1,zone=="I")
ps.pm.rad.lk1i <- prune_taxa(taxa_sums(ps.pm.rad.lk1i)!=0,ps.pm.rad.lk1i)
ps.pm.rad.lk1o <- subset_samples(ps.pm.rad.lk1,zone=="O")
ps.pm.rad.lk1o <- prune_taxa(taxa_sums(ps.pm.rad.lk1o)!=0,ps.pm.rad.lk1o)
ps.pm.rad.uk1i <- subset_samples(ps.pm.rad.uk1,zone=="I")
ps.pm.rad.uk1i <- prune_taxa(taxa_sums(ps.pm.rad.uk1i)!=0,ps.pm.rad.uk1i)
ps.pm.rad.uk1o <- subset_samples(ps.pm.rad.uk1,zone=="O")
ps.pm.rad.uk1o <- prune_taxa(taxa_sums(ps.pm.rad.uk1o)!=0,ps.pm.rad.uk1o)
ps.pm.rad.uk2i <- subset_samples(ps.pm.rad.uk2,zone=="I")
ps.pm.rad.uk2i <- prune_taxa(taxa_sums(ps.pm.rad.uk2i)!=0,ps.pm.rad.uk2i)
ps.pm.rad.uk2o <- subset_samples(ps.pm.rad.uk2,zone=="O")
ps.pm.rad.uk2o <- prune_taxa(taxa_sums(ps.pm.rad.uk2o)!=0,ps.pm.rad.uk2o)3.4.2 Time - LK1 rads
##time - inshore
seq.rad.lk1i <- data.frame(ps.pm.rad.lk1i@otu_table)
samdf.rad.lk1i <- data.frame(ps.pm.rad.lk1i@sam_data)
dist.rad.lk1i <- vegdist(seq.rad.lk1i)
bet.rad.lk1i <- betadisper(dist.rad.lk1i,samdf.rad.lk1i$time)
#anova(bet.ps.pm)
permutest(bet.rad.lk1i,pairwise=TRUE,permutations=999) ##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.90197 0.45098 31.687 999 0.001 ***
## Residuals 21 0.29888 0.01423
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 1.0000e-03 0.680
## Post 1.6560e-05 0.001
## Pre 6.7417e-01 4.1092e-05##Post different than other two
# Response: Distances
# Df Sum Sq Mean Sq F N.Perm Pr(>F)
# Groups 2 0.90197 0.45098 31.687 999 0.001 ***
# Residuals 21 0.29888 0.01423
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Pairwise comparisons:
# (Observed p-value below diagonal, permuted p-value above diagonal)
# Irma Post Pre
# Irma 1.0000e-03 0.664
# Post 1.6560e-05 0.001
# Pre 6.7417e-01 4.1092e-05
plot(bet.rad.lk1i) 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.rad.lk1i ~ time, data = samdf.rad.lk1i, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 1.1086 0.2633 3.7528 0.007 **
## Residual 21 3.1019 0.7367
## Total 23 4.2105 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Pre vs Irma 2.630912 0.1798187 0.024 0.024
## 2 Pre vs Post 3.634913 0.1950593 0.033 0.033
## 3 Irma vs Post 4.034110 0.2119411 0.013 0.013##all different
# pairs F.Model R2 p.value p.adjusted
# 1 Pre vs Irma 2.630912 0.1798187 0.036 0.036
# 2 Pre vs Post 3.634913 0.1950593 0.033 0.033
# 3 Irma vs Post 4.034110 0.2119411 0.024 0.024
##time - offshore
seq.rad.lk1o <- data.frame(ps.pm.rad.lk1o@otu_table)
samdf.rad.lk1o <- data.frame(ps.pm.rad.lk1o@sam_data)
dist.rad.lk1o <- vegdist(seq.rad.lk1o)
bet.rad.lk1o <- betadisper(dist.rad.lk1o,samdf.rad.lk1o$time)
#anova(bet.ps.pm)
permutest(bet.rad.lk1o,pairwise=TRUE,permutations=999)#ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.05765 0.028824 0.3725 999 0.674
## Residuals 18 1.39299 0.077388
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.35000 0.779
## Post 0.34828 0.606
## Pre 0.77610 0.61360
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.rad.lk1o ~ time, data = samdf.rad.lk1o, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 0.5146 0.09471 0.9416 0.453
## Residual 18 4.9186 0.90529
## Total 20 5.4332 1.00000## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Irma vs Pre 0.7328966 0.05755930 0.355 0.355
## 2 Irma vs Post 1.0368385 0.07953144 0.352 0.352
## 3 Pre vs Post 1.0230342 0.07855575 0.303 0.3033.4.3 Time - UK1 rads
##time
seq.rad.uk1i <- data.frame(ps.pm.rad.uk1i@otu_table)
samdf.rad.uk1i <- data.frame(ps.pm.rad.uk1i@sam_data)
dist.rad.uk1i <- vegdist(seq.rad.uk1i)
bet.rad.uk1i <- betadisper(dist.rad.uk1i,samdf.rad.uk1i$time)
#anova(bet.ps.pm)
permutest(bet.rad.uk1i,pairwise=TRUE,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.009434 0.0094338 1.1017 999 0.303
## Residuals 10 0.085631 0.0085631
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Pre
## Irma 0.306
## Pre 0.3186
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.rad.uk1i ~ time, data = samdf.rad.uk1i, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 1 0.12448 0.25484 3.4199 0.026 *
## Residual 10 0.36399 0.74516
## Total 11 0.48847 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Pre vs Irma 3.41993 0.2548396 0.033 0.0333.4.4 Time - UK2 rads
##time - inshore
seq.rad.uk2i <- data.frame(ps.pm.rad.uk2i@otu_table)
samdf.rad.uk2i <- data.frame(ps.pm.rad.uk2i@sam_data)
dist.rad.uk2i <- vegdist(seq.rad.uk2i)
bet.rad.uk2i <- betadisper(dist.rad.uk2i,samdf.rad.uk2i$time)
#anova(bet.ps.pm)
permutest(bet.rad.uk2i,pairwise=TRUE,permutations=999)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.22544 0.112720 5.305 999 0.015 *
## Residuals 21 0.44621 0.021248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.039000 0.225
## Post 0.033448 0.018
## Pre 0.222245 0.023323##Pre & irma are different
# Irma Post Pre
# Irma 0.026000 0.218
# Post 0.033448 0.016
# Pre 0.222245 0.023323
plot(bet.rad.uk2i) #post less constrained
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.rad.uk2i ~ time, data = samdf.rad.uk2i, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 0.8178 0.21125 2.8123 0.001 ***
## Residual 21 3.0534 0.78875
## Total 23 3.8712 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Pre vs Irma 1.514262 0.09169421 0.150 0.150
## 2 Pre vs Post 2.938851 0.19672535 0.008 0.008
## 3 Irma vs Post 3.344836 0.18233122 0.002 0.002##post different than other two
# pairs F.Model R2 p.value p.adjusted
# 1 Pre vs Irma 1.514262 0.09169421 0.161 0.161
# 2 Pre vs Post 2.938851 0.19672535 0.004 0.004
# 3 Irma vs Post 3.344836 0.18233122 0.005 0.005
##time - offshore
seq.rad.uk2o <- data.frame(ps.pm.rad.uk2o@otu_table)
samdf.rad.uk2o <- data.frame(ps.pm.rad.uk2o@sam_data)
dist.rad.uk2o <- vegdist(seq.rad.uk2o)
bet.rad.uk2o <- betadisper(dist.rad.uk2o,samdf.rad.uk2o$time)
#anova(bet.ps.pm)
permutest(bet.rad.uk2o,pairwise=TRUE,permutations=999)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.03579 0.017893 0.2306 999 0.805
## Residuals 23 1.78482 0.077601
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Irma Post Pre
## Irma 0.65000 0.562
## Post 0.66724 0.750
## Pre 0.56121 0.74220##Pre & irma are different
# Irma Post Pre
# Irma 0.026000 0.218
# Post 0.033448 0.016
# Pre 0.222245 0.023323
plot(bet.rad.uk2o) #post less constrained
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = seq.rad.uk2o ~ time, data = samdf.rad.uk2o, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## time 2 0.7432 0.10812 1.3941 0.178
## Residual 23 6.1308 0.89188
## Total 25 6.8740 1.00000## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead## pairs F.Model R2 p.value p.adjusted
## 1 Post vs Pre 2.0874958 0.10936457 0.081 0.081
## 2 Post vs Irma 1.6509889 0.10548783 0.118 0.118
## 3 Pre vs Irma 0.4521508 0.02926135 0.682 0.6823.5 Plot by time - rads
3.5.1 PCOA LK1 rads
ps.pm.rad.lk1i.no0 <- prune_taxa(taxa_sums(ps.pm.rad.lk1i)!=0,ps.pm.rad.lk1i)
ord.pm.rad.lk1i.no0 <- ordinate(ps.pm.rad.lk1i.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.rad.lk1i.no0, ord.pm.rad.lk1i.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()
ps.pm.rad.lk1o.no0 <- prune_taxa(taxa_sums(ps.pm.rad.lk1o)!=0,ps.pm.rad.lk1o)
ord.pm.rad.lk1o.no0 <- ordinate(ps.pm.rad.lk1o.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.rad.lk1o.no0, ord.pm.rad.lk1o.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
3.5.2 PCOA UK1 rads
ps.pm.rad.uk1i.no0 <- prune_taxa(taxa_sums(ps.pm.rad.uk1i)!=0,ps.pm.rad.uk1i)
ord.pm.rad.uk1i.no0 <- ordinate(ps.pm.rad.uk1i.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.rad.uk1i.no0, ord.pm.rad.uk1i.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()
ps.pm.rad.uk1o.no0 <- prune_taxa(taxa_sums(ps.pm.rad.uk1o)!=0,ps.pm.rad.uk1o)
ord.pm.rad.uk1o.no0 <- ordinate(ps.pm.rad.uk1o.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.rad.uk1o.no0, ord.pm.rad.uk1o.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()
3.5.3 PCOA UK2 rads
ps.pm.rad.uk2i.no0 <- prune_taxa(taxa_sums(ps.pm.rad.uk2i)!=0,ps.pm.rad.uk2i)
ord.pm.rad.uk2i.no0 <- ordinate(ps.pm.rad.uk2i.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.rad.uk2i.no0, ord.pm.rad.uk2i.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()
ps.pm.rad.uk2o.no0 <- prune_taxa(taxa_sums(ps.pm.rad.uk2o)!=0,ps.pm.rad.uk2o)
ord.pm.rad.uk2o.no0 <- ordinate(ps.pm.rad.uk2o.no0,"PCoA",distance="bray")
plot_ordination(ps.pm.rad.uk2o.no0, ord.pm.rad.uk2o.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
#facet_grid(site~zone)+
theme_cowplot()## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
4 Script graveyard
Type profile business
setwd("/Volumes/GoogleDrive/My Drive/Flirma/flirma_networks/fl_its2")
#ps <- readRDS("ps.its2.rds")
#ps.rel <- transform_sample_counts(ps, function(x) x / sum(x))
#saveRDS(ps.rel,"ps.its2.rel.rds")
ps1 <- readRDS("ps.its2.rel.rds")
ps <- subset_samples(ps1,site!="UK3")
ps.sid <- subset_samples(ps,spp=="ssid")
ps.sid.pre <- subset_samples(ps.sid,time=="Pre")
ps.sid.irm <- subset_samples(ps.sid,time=="Irma")
ps.sid.pos <- subset_samples(ps.sid,time=="Post")
ps.sid.pre.lk1 <- subset_samples(ps.sid.pre,site=="LK1")
ps.sid.pre.uk1 <- subset_samples(ps.sid.pre,site=="UK1")
ps.sid.pre.uk2 <- subset_samples(ps.sid.pre,site=="UK2")
ps.sid.irm.lk1 <- subset_samples(ps.sid.irm,site=="LK1")
ps.sid.irm.uk1 <- subset_samples(ps.sid.irm,site=="UK1")
ps.sid.irm.uk2 <- subset_samples(ps.sid.irm,site=="UK2")
ps.sid.pos.lk1 <- subset_samples(ps.sid.pos,site=="LK1")
ps.sid.pos.uk1 <- subset_samples(ps.sid.pos,site=="UK1")
ps.sid.pos.uk2 <- subset_samples(ps.sid.pos,site=="UK2")
##for time stats
ps.sid.lk1 <- subset_samples(ps.sid,site=="LK1")
ps.sid.uk1 <- subset_samples(ps.sid,site=="UK1")
ps.sid.uk2 <- subset_samples(ps.sid,site=="UK2")
ps.sid.lk1i <- subset_samples(ps.sid.lk1,zone=="I")
ps.sid.lk1o <- subset_samples(ps.sid.lk1,zone=="O")
ps.sid.uk1i <- subset_samples(ps.sid.uk1,zone=="I")
ps.sid.uk1o <- subset_samples(ps.sid.uk1,zone=="O")
ps.sid.uk2i <- subset_samples(ps.sid.uk2,zone=="I")
ps.sid.uk2o <- subset_samples(ps.sid.uk2,zone=="O")4.0.1 Plots by time
ps.sid.lk1i.no0 <- prune_taxa(taxa_sums(ps.sid.lk1i)!=0,ps.sid.lk1i)
ord.sid.lk1i.no0 <- ordinate(ps.sid.lk1i.no0,"PCoA",distance="bray")
plot_ordination(ps.sid.lk1i.no0, ord.sid.lk1i.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
theme_cowplot()
ps.sid.lk1o.no0 <- prune_taxa(taxa_sums(ps.sid.lk1o)!=0,ps.sid.lk1o)
ord.sid.lk1o.no0 <- ordinate(ps.sid.lk1o.no0,"PCoA",distance="bray")
plot_ordination(ps.sid.lk1o.no0, ord.sid.lk1o.no0, color ="time")+
geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
theme_cowplot()4.0.2 Reef zones
ord.sid.pre.lk1 <- ordinate(ps.sid.pre.lk1,"PCoA",distance="bray")
plot_ordination(ps.sid.pre.lk1, ord.sid.pre.lk1, color ="zone")+
#geom_point(alpha=0.5)+
#scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
#scale_shape_manual(values=c(16,15))+
stat_ellipse()+
theme_cowplot()
#facet_grid(site~zone)
#ggsave(filename="pcoa.types.site.pdf")5 Stats - type profiles
5.1 Setup
5.2 Subset
ps.sid <- subset_samples(ps,spp="ssid")
ps.sid <- transform_sample_counts(ps.sid, function(x) x / sum(x))
##for in v off stats
ps.sid.pre <- subset_samples(ps.sid,time=="Pre")
ps.sid.irm <- subset_samples(ps.sid,time=="Irma")
ps.sid.pos <- subset_samples(ps.sid,time=="Post")
ps.sid.pre.lk1 <- subset_samples(ps.sid.pre,site=="LK1")
ps.sid.pre.uk1 <- subset_samples(ps.sid.pre,site=="UK1")
ps.sid.pre.uk2 <- subset_samples(ps.sid.pre,site=="UK2")
ps.sid.irm.lk1 <- subset_samples(ps.sid.irm,site=="LK1")
ps.sid.irm.uk1 <- subset_samples(ps.sid.irm,site=="UK1")
ps.sid.irm.uk2 <- subset_samples(ps.sid.irm,site=="UK2")
ps.sid.pos.lk1 <- subset_samples(ps.sid.pos,site=="LK1")
ps.sid.pos.uk1 <- subset_samples(ps.sid.pos,site=="UK1")
ps.sid.pos.uk2 <- subset_samples(ps.sid.pos,site=="UK2")
##for time stats
ps.sid.lk1 <- subset_samples(ps.sid,site=="LK1")
ps.sid.uk1 <- subset_samples(ps.sid,site=="UK1")
ps.sid.uk2 <- subset_samples(ps.sid,site=="UK2")
ps.sid.lk1i <- subset_samples(ps.sid.lk1,zone=="I")
ps.sid.lk1o <- subset_samples(ps.sid.lk1,zone=="O")
ps.sid.uk1i <- subset_samples(ps.sid.uk1,zone=="I")
ps.sid.uk1o <- subset_samples(ps.sid.uk1,zone=="O")
ps.sid.uk2i <- subset_samples(ps.sid.uk2,zone=="I")
ps.sid.uk2o <- subset_samples(ps.sid.uk2,zone=="O")5.3 Adonisii - sids
5.3.1 Reef zones - LK1
##reef zone pre
seq.sid.pre.lk1 <- data.frame(ps.sid.pre.lk1@otu_table)
samdf.sid.pre.lk1 <- data.frame(ps.sid.pre.lk1@sam_data)
dist.sid.pre.lk1 <- vegdist(seq.sid.pre.lk1)
##zone
bet.sid.pre.lk1 <- betadisper(dist.sid.pre.lk1,samdf.sid.pre.lk1$zone)
#anova(bet.ps)
permutest(bet.sid.pre.lk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.pre.lk1) #ns
adonis2(seq.sid.pre.lk1 ~ zone, data=samdf.sid.pre.lk1, permutations=999) #ns
##reef zone irm
seq.sid.irm.lk1 <- data.frame(ps.sid.irm.lk1@otu_table)
samdf.sid.irm.lk1 <- data.frame(ps.sid.irm.lk1@sam_data)
dist.sid.irm.lk1 <- vegdist(seq.sid.irm.lk1)
##zone
bet.sid.irm.lk1 <- betadisper(dist.sid.irm.lk1,samdf.sid.irm.lk1$zone)
#anova(bet.ps)
permutest(bet.sid.irm.lk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.irm.lk1) #ns
adonis2(seq.sid.irm.lk1 ~ zone, data=samdf.sid.irm.lk1, permutations=999)
#adonis2(formula = seq.sid.irm.lk1 ~ zone, data = samdf.sid.irm.lk1, permutations = 999)
# Df SumOfSqs R2 F Pr(>F)
# zone 1 0.8368 0.06787 1.966 0.038 *
# Residual 27 11.4929 0.93213
# Total 28 12.3298 1.00000
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##reef zone post
seq.sid.pos.lk1 <- data.frame(ps.sid.pos.lk1@otu_table)
samdf.sid.pos.lk1 <- data.frame(ps.sid.pos.lk1@sam_data)
dist.sid.pos.lk1 <- vegdist(seq.sid.pos.lk1)
##zone
bet.sid.pos.lk1 <- betadisper(dist.sid.pos.lk1,samdf.sid.pos.lk1$zone)
#anova(bet.ps)
permutest(bet.sid.pos.lk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.pos.lk1) #ns
adonis2(seq.sid.pos.lk1 ~ zone, data=samdf.sid.pos.lk1, permutations=999) #ns5.3.2 Reef zones - UK1
##reef zone pre
seq.sid.pre.uk1 <- data.frame(ps.sid.pre.uk1@otu_table)
samdf.sid.pre.uk1 <- data.frame(ps.sid.pre.uk1@sam_data)
dist.sid.pre.uk1 <- vegdist(seq.sid.pre.uk1)
##zone
bet.sid.pre.uk1 <- betadisper(dist.sid.pre.uk1,samdf.sid.pre.uk1$zone)
#anova(bet.ps)
permutest(bet.sid.pre.uk1,pairwise=FALSE,permutations=999)
plot(bet.sid.pre.uk1)
# Df SumOfSqs R2 F Pr(>F)
# zone 1 0.914 0.06604 2.0507 0.023 *
# Residual 29 12.925 0.93396
# Total 30 13.839 1.00000
adonis2(seq.sid.pre.uk1 ~ zone, data=samdf.sid.pre.uk1, permutations=999)
# adonis2(formula = seq.sid.pre.uk1 ~ zone, data = samdf.sid.pre.uk1, permutations = 999)
# Df SumOfSqs R2 F Pr(>F)
# zone 1 0.914 0.06604 2.0507 0.021 *
# Residual 29 12.925 0.93396
# Total 30 13.839 1.00000
# ---
##reef zone irm
seq.sid.irm.uk1 <- data.frame(ps.sid.irm.uk1@otu_table)
samdf.sid.irm.uk1 <- data.frame(ps.sid.irm.uk1@sam_data)
dist.sid.irm.uk1 <- vegdist(seq.sid.irm.uk1)
##zone
bet.sid.irm.uk1 <- betadisper(dist.sid.irm.uk1,samdf.sid.irm.uk1$zone)
#anova(bet.ps)
permutest(bet.sid.irm.uk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.irm.uk1) #ns
adonis2(seq.sid.irm.uk1 ~ zone, data=samdf.sid.irm.uk1, permutations=999)
# adonis2(formula = seq.sid.irm.uk1 ~ zone, data = samdf.sid.irm.uk1, permutations = 999)
# Df SumOfSqs R2 F Pr(>F)
# zone 1 1.3958 0.09887 3.2915 0.002 **
# Residual 30 12.7216 0.90113
# Total 31 14.1174 1.00000
##reef zone post
seq.sid.pos.uk1 <- data.frame(ps.sid.pos.uk1@otu_table)
samdf.sid.pos.uk1 <- data.frame(ps.sid.pos.uk1@sam_data)
dist.sid.pos.uk1 <- vegdist(seq.sid.pos.uk1)
##zone
bet.sid.pos.uk1 <- betadisper(dist.sid.pos.uk1,samdf.sid.pos.uk1$zone)
#anova(bet.ps)
permutest(bet.sid.pos.uk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.pos.uk1) #ns
adonis2(seq.sid.pos.uk1 ~ zone, data=samdf.sid.pos.uk1, permutations=999) #ns5.3.3 Reef zone - UK2
5.3.4 Time - LK1 sids
##time - inshore
seq.sid.lk1i <- data.frame(ps.sid.lk1i@otu_table)
samdf.sid.lk1i <- data.frame(ps.sid.lk1i@sam_data)
dist.sid.lk1i <- vegdist(seq.sid.lk1i)
bet.sid.lk1i <- betadisper(dist.sid.lk1i,samdf.sid.lk1i$time)
#anova(bet.ps)
permutest(bet.sid.lk1i,pairwise=TRUE,permutations=999)
#post different than other two, looks more constrained
# Permutation test for homogeneity of multivariate dispersions
# Permutation: free
# Number of permutations: 999
#
# Response: Distances
# Df Sum Sq Mean Sq F N.Perm Pr(>F)
# Groups 2 0.05045 0.0252251 4.4942 999 0.016 *
# Residuals 46 0.25819 0.0056128
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Pairwise comparisons:
# (Observed p-value below diagonal, permuted p-value above diagonal)
# Irma Post Pre
# Irma 0.0240000 0.349
# Post 0.0239826 0.006
# Pre 0.3359886 0.0082901
plot(bet.sid.lk1i)
adonis2(seq.sid.lk1i ~ time, data=samdf.sid.lk1i, permutations=999) #p<0.001***
pairwise.adonis(seq.sid.lk1i, factors=samdf.sid.lk1i$time, permutations=999)
##pre diff than post
# pairs F.Model R2 p.value p.adjusted
# 1 Pre vs Irma 1.273944 0.04208053 0.272 0.272
# 2 Pre vs Post 3.211953 0.08869871 0.001 0.001
# 3 Irma vs Post 1.507718 0.04785234 0.124 0.124
##time - offshore
seq.sid.lk1o <- data.frame(ps.sid.lk1o@otu_table)
samdf.sid.lk1o <- data.frame(ps.sid.lk1o@sam_data)
dist.sid.lk1o <- vegdist(seq.sid.lk1o)
bet.sid.lk1o <- betadisper(dist.sid.lk1o,samdf.sid.lk1o$time)
#anova(bet.ps)
permutest(bet.sid.lk1o,pairwise=TRUE,permutations=999)#ns
plot(bet.sid.lk1o)
adonis2(seq.sid.lk1o ~ time, data=samdf.sid.lk1o, permutations=999) #p<0.01**
pairwise.adonis(seq.sid.lk1o, factors=samdf.sid.lk1o$time, permutations=999)
#post different than other two
# pairs F.Model R2 p.value p.adjusted
# 1 Irma vs Pre 1.549658 0.05072586 0.132 0.132
# 2 Irma vs Post 2.674018 0.08442309 0.001 0.001
# 3 Pre vs Post 2.582747 0.07926733 0.014 0.0145.3.5 Time - UK1 sids
##time
seq.sid.uk1i <- data.frame(ps.sid.uk1i@otu_table)
samdf.sid.uk1i <- data.frame(ps.sid.uk1i@sam_data)
dist.sid.uk1i <- vegdist(seq.sid.uk1i)
bet.sid.uk1i <- betadisper(dist.sid.uk1i,samdf.sid.uk1i$time)
#anova(bet.ps)
permutest(bet.sid.uk1i,pairwise=TRUE,permutations=999) #ns
plot(bet.sid.uk1i)
adonis2(seq.sid.uk1i ~ time, data=samdf.sid.uk1i, permutations=999) #p<0.01**
pairwise.adonis(seq.sid.uk1i, factors=samdf.sid.uk1i$time, permutations=999)
# pairs F.Model R2 p.value p.adjusted
# 1 Post vs Pre 1.903144 0.07347155 0.026 0.026
# 2 Post vs Irma 2.306884 0.08149552 0.022 0.022
# 3 Pre vs Irma 1.898341 0.05600101 0.047 0.047
#all different5.3.6 Time - UK2 sids
##time
seq.sid.uk2i <- data.frame(ps.sid.uk2i@otu_table)
samdf.sid.uk2i <- data.frame(ps.sid.uk2i@sam_data)
dist.sid.uk2i <- vegdist(seq.sid.uk2i)
bet.sid.uk2i <- betadisper(dist.sid.uk2i,samdf.sid.uk2i$time)
#anova(bet.ps)
permutest(bet.sid.uk2i,pairwise=TRUE,permutations=999)
##Pre & irma are different
# Irma Post Pre
# Irma 0.439000 0.016
# Post 0.431974 0.187
# Pre 0.017649 0.188058
plot(bet.sid.uk2i) #pre less constrained
adonis2(seq.sid.uk2i ~ time, data=samdf.sid.uk2i, permutations=999) #p<0.05*
pairwise.adonis(seq.sid.uk2i, factors=samdf.sid.uk2i$time, permutations=999)
##post different than other two
# pairs F.Model R2 p.value p.adjusted
# 1 Pre vs Post 1.958299 0.05766776 0.028 0.028
# 2 Pre vs Irma 1.398358 0.03461423 0.164 0.164
# 3 Post vs Irma 1.613048 0.03971747 0.050 0.050ps.sid.rel <- transform_sample_counts(ps.sid, function(x) x / sum(x))
## Presence/absence
sppAbun <- as.data.frame(ps.sid.rel@otu_table)
sppAbun[sppAbun > 0] <- 1 #converts from abundance to P/A
sppAbun
rowSums(sppAbun)
morethan1 <- sppAbun[rowSums(sppAbun)>1,]
find.top.taxa <- function(x,taxa){
require(phyloseq)
top.taxa <- tax_glom(x, taxa)
otu <- otu_table(t(top.taxa)) # remove the transformation if using a merge_sample object
tax <- tax_table(top.taxa)
j<-apply(otu,1,which.max)
k <- j[!duplicated(j)]
l <- data.frame(tax[k,])
m <- data.frame(otu[,k])
s <- as.name(taxa)
colnames(m) = l[,taxa]
n <- colnames(m)[apply(m,1,which.max)]
m[,taxa] <- n
return(m)
}
find.top.taxa(top.pdy,"Phylum")
#chi-square, top vs side
chisq.test(SymPortalDataK_1$variable, SymPortalDataK_1$label, correct=FALSE)
#Pearson's Chi-squared test
#data: SymPortalDataK_1$variable and SymPortalDataK_1$label
#X-squared = 5.5824, df = 6, p-value = 0.47165.4 Find “most frequent” sym
Just ran once to get the info
find.top.asv <- function(x,num){
require(phyloseq)
require(magrittr)
otu <- otu_table(x)
tax <- tax_table(x)
j1 <- apply(otu,1,sort,index.return=T, decreasing=T) # modifying which.max to return a list of sorted index
j2 <- lapply(j1,'[[',"ix") # select for index
l <- data.frame(unique(tax@.Data[unlist(j2),]))
m <- data.frame(otu@.Data[,unique(unlist(j2))])
n <- apply(m,1,sort,index.return=T, decreasing=T) %>%
lapply('[[',"ix") %>% # Extract index
lapply(head,n=num) # This to returns the top x tax
p <- list()
for(i in 1:length(n)){
p[[i]]<- colnames(m)[n[[i]]]
}
m$taxa <- p
return(m)
}
top.asvs <- find.top.asv(ps.sid,5)
top.asvs$taxa <- as.character(top.asvs$taxa)
#write.csv(top.asvs,file="top_sym.csv",row.names=TRUE)6 Flirma scripts scraps
6.1 Collapsing type profiles via phylogeny tree
Started with file: “188_20211214_03_DBV_20211215T030311.seqs.fasta” in the post-med seqs folder. Also started with the taxa info ‘Majority_ITS2_sequence", removed the slashes between types, & removed duplicates. That left 24 unique ITS2 types, put these in ’names.txt’. Want to extract the relevant sequences from the fasta file now
for word in $(cat names.txt); do echo $word; grep -w -A 1 $word postmedseqs.fasta >> output.fasta; doneImportant notes: produced a fasta file with ‘–’ in between some desired entries, removed that manually, but I’m sure there’s a better way than that. It also gave me an unwanted one (A4.3), took that out manually. Easy to do with 24 types, but would be harder with more than that.. processed file was named output_types_cleaned.fasta
Then went to ngphylogeny.fr to plot.
More notes: - A is fine, the majority sequences are pretty disparate. - The Bs are a total mess, will plot below - Cs are fine - Ds are fine
#facet_grid(site~zone)
#
# #devtools::install_github("david-barnett/microViz@0.9.3") # check 0.9.3 is the latest release
# library("microViz")
#
# ps.sid.lk1i.no0 %>%
# ps_mutate(
# time = as.numeric(time == "time")
# ) %>%
# ord_calc(
# constraints = c("time"),
# # method = "RDA", # Note: you can specify RDA explicitly, and it is good practice to do so, but microViz can guess automatically that you want an RDA here (helpful if you don't remember the name?)
# scale_cc = FALSE # doesn't make a difference
# ) %>%
# plot_ord(
# colour = "time", size = 2, alpha = 0.5, shape = "active",
# plot_taxa = 1:8
# )